
local lhs cog
local lowp p25
local highp p75
local med p50

capture program drop ccn_beta_delta_het_newsym
local subsample = "$subsample"

clear
local count=0
local controls="${CONTROLlist_`lhs'}"
local beta_controls "${BETACONTROLlist_hrs_cog_linint}"
local delta_controls "${DELTACONTROLlist_hrs_cog_linint}"
foreach rhs in $RHSlist { 
	foreach robust in $robust_list {
		forvalues K=$K_list {
			noisily display "LHS=`lhs', RHS=`rhs', SE=`robust', K=`K', $S_DATE, $S_TIME"
			use "$datapath/003_cluster_`subsample'_`lhs'.dta", clear
			local control_length = 0
			foreach control of varlist `controls' {
				local control_length = `control_length' + 1
			}		
			global control_length = `control_length'
			qui tab X_`K', generate(dumX_`K'_)

			
			*loop over types
			local model_types symmetric
			foreach model of local model_types {
				di "`model'"
				if ("`model'" == "het_tobit") {	
					local title_name = "Het. Tobit"
				}

				if ("`model'" == "het_uniform") {	
					local title_name "Het. Uniform"
				}
				if ("`model'" == "symmetric") {	
					local title_name "Symmetric"
				}

				*********************************************************************************************************
				*Look at all hours for a fixed set of skills
				*Uniform
				qui ccn_beta_delta_het_newsym `model' `K' `rhs' `lhs' `robust' "`controls'" het_tobit 1 "`beta_controls'" "`delta_controls'" 0
				local b0 = COEFFS[1,3] // constant
				local bhours = COEFFS[1,4] // hours
				local bafqt = COEFFS[1,5] // afqt
				local bint = COEFFS[1,6] // hrsafqt
				*noisily di "b0 is `b0', bhours is `bhours', bafqt is `bafqt', bint is `bint'"
				
				xtile afqt_mom_buck = afqt_mom, nq(4)
				forvalues ii=1(1)4 {
					qui sum afqt_mom if afqt_mom_buck == `ii'
					local mom_skill_`ii' = r(mean)
					qui sum `lhs' if afqt_mom_buck == `ii'
					local skill_`ii' = r(mean)
					qui sum hours if afqt_mom_buck == `ii' // & hours != 0 //condition on non-0 or no?
					local hours_`ii' = r(mean)
					local actual_beta_`ii' = (`b0'+ (`bhours'*`hours_`ii'') + (`bafqt'*`mom_skill_`ii'') + (`bint'*`mom_skill_`ii''*`hours_`ii''))
					local actual_effect_`ii' = (`b0'+ (`bhours'*`hours_`ii'') + (`bafqt'*`mom_skill_`ii'') + (`bint'*`mom_skill_`ii''*`hours_`ii'')) * `hours_`ii''
					noisily di "actual effect = `actual_effect_`ii'', actual_beta_`ii' = `actual_beta_`ii''"
					noisily di "mom_skill_`ii' = `mom_skill_`ii'', hours_`ii'= `hours_`ii'', skill_`ii'=`skill_`ii''"
				}
				
				*******
				rename hours hours_original
				gen hours = 0 in 1
				replace hours = hours[_n-1] + 100 in 2/26
				
				forvalues ii = 1(1)4 {
					gen pred_effect_`ii' = (`b0'+`bhours'*hours + `bafqt'*`mom_skill_`ii''+`bint'*`mom_skill_`ii''*hours)*hours - `actual_effect_`ii''
					gen pred_`ii' = `skill_`ii''+pred_effect_`ii'
				}

				gen pred_gap = pred_4 - pred_1
				gen pred_gap_change = (pred_effect_4 - pred_effect_1)

				gen zero_hour = (hours_original==0)
				sum zero_hour if afqt_mom_buck==1
				local value_label_1  = round(100*`r(mean)', 0.1)
				gen zero_hour_1 = 18 in 1

				label define zero_hour_1 18 "`value_label_1'", replace
				label values zero_hour_1 zero_hour_1

				sum zero_hour if afqt_mom_buck == 4
				local value_label_4  = round(100*`r(mean)', 0.1)
				gen zero_hour_4 = 15 in 1

				label define zero_hour_4 15 "`value_label_4'", replace
				label values zero_hour_4 zero_hour_4

					
				*High Low Skills with no gap, no 0 line, no histogram, Markers for the observed hours/skills
				forvalues ii=1(1)4 {
					gen hours_`ii' = .
					replace hours_`ii' = `hours_`ii'' if _n == 1
					gen skill_`ii' = .
					replace skill_`ii' = `skill_`ii'' if _n == 1
				}
				
				
				*Now Show All 4
				twoway ///
					(line pred_1 hours, sort lcolor(black)) ///
					(line pred_2 hours, sort lcolor(black) lpattern(dot)) ///
					(line pred_3 hours, sort lcolor(black)lpattern(shortdash_dot)) ///
					(line pred_4 hours, sort lpattern(dash) lcolor(black)) ///
					(scatter skill_1 hours_1, msymbol(circle_hollow) mcolor(black) msize(large) ) ///
					(scatter skill_2 hours_2, msymbol(circle_hollow) mcolor(black) msize(large) ) ///
					(scatter skill_3 hours_3, msymbol(circle_hollow) mcolor(black) msize(large) ) ///
					(scatter skill_4 hours_4, msymbol(circle_hollow) mcolor(black) msize(large)), scheme(s1mono) ///
					legend(order(1 "Quartile 1" 2 "Quartile 2" 3 "Quartile 3" 4 "Quartile 4") ///
					cols(2) region(lcolor(none))) ///
					ytitle("Skills (sd)") xtitle(Hours)
					graph export  "$figpath/total_skills_high_low_beta_delta_het_`lhs'_linint_`model'_all_quartiles.pdf", replace
	
					
			}

		}

	}
}





